PARAXIAL COUPLING OF PROPAGATING MODES IN THREE-DIMENSIONAL 
WAVEGUIDES WITH RANDOM BOUNDARIES 
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Abstract. We analyze long range wave propagation in three-dimensional random waveguides. The waves are trapped by 
top and bottom boundaries, but the medium is unbounded in the two remaining directions. We consider scalar waves, and 
motivated by applications in underwater acoustics, we take a pressure release boundary condition at the top surface and a rigid 
bottom boundary. The wave speed in the waveguide is known and smooth, but the top boundary has small random fluctuations 
that cause significant cumulative scattering of the waves over long distances of propagation. To quantify the scattering effects, 
we study the evolution of the random amplitudes of the waveguide modes. We obtain that in the long range limit they satisfy 
a system of paraxial equations driven by a Brownian field. We use this system to estimate three important mode-dependent 
scales: the scattering mean free path, the cross-range decoherence length and the decoherence frequency. Understanding these 
scales is important in imaging and communication problems, because they encode the cumulative scattering effects in the wave 
field measured by remote sensors. As an application of the theory, we analyze time reversal and coherent interfero metric imaging 
in strong cumulative scattering regimes. 
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1. Introduction. We study long range scalar (acoustic) wave propagation in a three-dimensional 
waveguide. The setup is illustrated in Figure and it is motivated by problems in underwater acous- 
Ph tics. We denote by z € R the range, the main direction of propagation of the waves. The medium is 

unbounded in the cross-range direction x £ R, but it is confined in depth y by two boundaries which trap 
f-H the waves, thus creating the waveguide effect. 

The acoustic pressure field is denoted by pit, x, y, z), and it satisfies the wave equation 
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V c 2( y ) 



p(t,x,y,z) = f(t,x,y)d(z), y € [0, T(x, z)], i,z6l, t > 0, (1.1) 



t-H in a medium with wave speed c(y). The excitation is due to a source located in the plane z = 0, emitting 

the pulse f(t). The medium is quiescent before the source excitation, 

p(t,x,y,z) = 0, <<0. (1.2) 

The bottom of the waveguide is assumed rigid 

d y p(t,x,y = 0,z)=0, (1.3) 

C^l and we take a pressure release boundary condition at the perturbed top boundary 

p(t,x,y = T(x,z),z) = 0. (1.4) 

Perturbed means that the boundary y = T(x, z) has small fluctuations around the mean depth T>, 

^ \T(x,z)-V\<ZV. (1.5) 

We choose this setup for simplicity. The results extend readily to other boundary conditions and to fluctuat- 
ing bottoms. Such boundaries were considered recently in [Tllll). in two-dimensional waveguides. Extensions 
to media with small (x, y, z)-dependent fluctuations of the wave speed can also be made using the techniques 
developed in [TJ [I M EH • 

The goal of our study is to quantify the effect of scattering at the surface. Because in applications it 
is not feasible to know the boundary fluctuations in detail, we model them with a random process. The 



solution p(t, x, y, z) of equations ( 1.1 1-( 1.4 1 is therefore a random field, and we describe in detail its statistics 
at long ranges, where cumulative scattering is significant. We use the results for two applications: time 
reversal and sensor array imaging. 
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y = V 



y = Q 

Fig. 1.1. Schematic of the problem setup. The system of coordinates has range origin z = at the source. The rigid 
bottom boundary y = is assumed flat and the pressure release top boundary has fluctuations around the value y = T>. The 
cross-range x and the range z are unbounded, that is (x,z) £ R 2 . 



Our method of solution uses a change of coordinates to straighten the boundary. The transformed 
problem has a simple geometry but a randomly perturbed differential operator. Its solution is given by a 
superposition of propagating and evanescent waveguide modes, with random amplitudes. We show that in 
the long range limit these amplitudes satisfy a system of paraxial equations that are driven by a Brownian 
field. The detailed characterization of the statistics of p(t,x,y,z) follows from this system. It involves the 
calculation of the mode-dependent scattering mean free path, which is the distance over which the modes 
lose coherence; the mode-dependent decoherence length, which is the cross-range offset over which the mode 
amplitudes decorrelate; and the mode-dependent decoherence frequency, which is the frequency offset over 
which the mode amplitudes decorrelate. These scales are important in studies of time reversal and imaging, 
because they dictate the resolution of focusing and the robustness (statistical stability) of the results with 
respect to realizations of the random fluctuations of the boundary. 

The paper is organized as follows: We begin in section [2] with the description of the reference pressure 
field p (t, x, y, z) in ideal waveguides with planar boundaries. The random field p{t, x, y, z) derived in section 
[3] may be viewed as a perturbation of p (t, x, y, z), in the sense that it is decomposed in the same waveguide 
modes. However, the amplitudes of the modes are random and coupled. Because the fluctuations of the 
boundary are small, we consider in section [4] a long range limit, so that we can observe significant cumulative 
scattering. The statistics of the wave field at such long ranges is described in section [5] The results are 
summarized in section [6] and are used in sections [7] and [H] for analyzing time reversal and imaging with sensor 
arrays. We end with a summary in section [9] 

2. Wave propagation in ideal waveguides. The pressure field in ideal waveguides, with planar 
boundaries, is given by 

/oo 
— p (u),x,y,z)e~ iUt , (2.1) 

with Fourier coefficients satisfying a separable problem for the Hclmholtz equation 



5 2 , Q 2 , c,2 , w 
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p (w,x,y,z) = f{u>,x,y)5{z), \ u -uo\<j, {x 7 z)eR 2 , y e (0,2?), (2.2) 



with boundary conditions 

d y p (uj, x,y = 0,z)= p (uj, x,y = V,z) = 0, (2.3) 
and outgoing radiation conditions at \/x 2 + z 2 — > oo. The Fourier transform of the source 

f(u,x,y)= f°° dtf(t,x,y)e iut , (2.4) 

J — oo 

is compactly supported in [ujq — B/2, ujq + B/2], for any x and y. Here ujq is the central frequency and B is 
the bandwidth. 
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2.1. Propagating and evanescent modes. The solution of the Hclmholtz equation (2.2 1 is a super- 
position of N(u>) propagating modes, and infinitely many evanescent ones, 



N(uj) oo 

Po(u,x,y,z) = 22 <f>j(u,y)u jt0 (u,x,z) + ^2 <t>Au,y)vj,a{u,x,z). 

3 = 1 j=N(u)+l 



(2.5) 



The decomposition is in the L 2 (0,T>) orthonormal basis of the eigenfunctions cf>j(uj,y) of the self-adjoint 
differential operator in y, 



d 2 + — 
v c 2 {y) 



t j (u,V) = B y <l> j (w ) 0)=0, j = 1,2, 



with eigenvalues Xj(to) that are simple [16j . 

To simplify the analysis, we assume in this paper that the wave speed is homogeneous 

c(y) = c . (2.6) 

The results for variable wave speeds are similar in all the essential aspects. The simplification brought by 
(2.6) amounts to having explicit expressions of the eigenfunctions, which are independent of the frequency 



i(y) 



cos 



1" [J 



2JV 



The eigenvalues are 



l\ 2 



j 



where k = co/c is the wavenumber, and only the first N(oj) of them are non-negative 

kV 1 

N(u) 



(2.7) 



(2.8) 



(2.9) 



The notation |_ J stands for the integer part. We suppose for simplicity that N{oS) remains constant in the 
bandwidth [wq — B/2,ujq + B/2], and write from now on N(u>) — N. 

The propagating components in (2.5) satisfy the two-dimensional Hclmholtz equation 



[dl + d 2 + $(w)] u j>0 (u, x, z) = Fj(u, x)S(z), j = l,...,N, 



(2.10) 



with outgoing, radiation conditions at \J x 2 + z 2 — > oo. The evanescent components solve 

[d 2 x + dl-p 2 3 {uj)]v h0 {uj,x,z)^F 3 {u J: x)5{z), j>N, (2.11) 

with decay condition Vj^ (uj, x, z) — > at \J x 2 + z 2 — > oo. Here we introduced the coefficients of the source 
profile in the basis of the eigenfunctions 



Fj(uj,x) 



dy<t>j(y)f(u,x,y), j>l, 



and the mode wavenumbers 



V 



kV\ 2 



-I -Ij- 



V 7T / 



1\ 2 



3> 1. 



(2.12) 



(2.13) 



We assume that none of the /3j(w) vanishes in the bandwidth, so that there are no standing waves. That is 
to say, 

kV 



N + a{u)-\, a(w)e(0,l) for all ui e [ui ~ B/2,w + B/2]. 
it 2 



(2.14) 



2.2. The paraxial regime. We now introduce the paraxial scaling for the ideal waveguide, with the 
source emitting a beam that propagates along the z-axis. As we show below, this happens when the cross- 
range profile of the source is larger than the wavelength. The source generates a quasi-plane wave, with 
slowly varying envelope satisfying a Schrodinger-like equation. 

Explicitly, we assume that the source is of the form 



f E (t,x,y) = f{t,ex,y) 



(2.15) 



where e is a small dimensionless parameter defined as the ratio of the central wavelength Ao and the transverse 
width ro of the source. Standard diffraction theory gives that the Rayleigh length for a beam with initial 
width 7" = A /e is of the order of 

ro/Ao - A /e 2 . 

The Rayleigh length is the distance along the z axis from the beam waist to the place where the beam area 
is doubled by diffraction. Therefore, we look at the wavefield at 0(e~ 1 ) cross-range scales, similar to ro, and 
at 0(s~ 2 ) range scale, similar to the Rayleigh length. We rename the field in this scaling as 



P £ o(t,X,y,Z) =p (t, —,y, 



The Fourier coefficients of (2.161 are given by the scaled version of (2.5) 



N(ui) 



3=1 j=N(u) + l 



(2.16) 



(2.17) 



with propagating mode amplitudes u e satisfying the scaled equation (2.10), with the source replaced by 
Fj(u>,ex — X). They can be written as 



m {u,x, z) 



dX 



x-x' z 



in terms of the outgoing Green's function 



Here Hq is the Hankel function of the first kind, and because e <C 1, we can use its asymptotic form for a 
scaled range Z > 



1 rrW 
4 H ° 



. (x - x 1 ) 2 z 2 



2? 



(x-x>) 2 



1/2 



exp 



. (X - X') 2 z 2 



exp < ifij{u)) 



2 y 2nPj(w)Z 
The propagating components of the wave field become 

u| j0 (w, X, Z) sa aj i0 (u;, X, Z) exp 

with 

1 / i 



(X - X') 2 
2Z 
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a jt0 (uj,X, Z) = - 



2 V 2n0j(w)Z 



dX' exp 



i^){x - x'y 

2Z 



Fj(uJ,X'), 



(2.18) 



for j = 1, ... , N. The evanescent components are obtained similarly from (2.11 1 



V s (uj,X,Z) w ej „(c«;, JC,Z)exp 



and 



ej t0 (u),X, Z) 



2irPj(tjj)Z 



dX' exp 



2Z 



(2.19) 



for j > JV + 1. These modes are exponentially damped and can be neglected. 
In summary, the paraxial approximation of the wave field is given by 



N 



4=1 



(2.20) 



It is a superposition of forward going modes with complex valued amplitudes a^ Q given by (2.18), and solving 
the paraxial equations 



[2ifa{w)dz + d 2 x ] a jt0 (u, X, Z) = 0, j = 1, 



with initial conditions 



a,j (u},X,Z = 0) 



F^X), j = l, 



(2.21) 



(2.22) 



3. Wave propagation in random waveguides. In this section we consider a waveguide with fluctu- 
ating boundary and analyze the wave field under the following scaling assumptions: 
1. The transverse width r of the source and the central wavelength X satisfy 



r =e 1 X , 



(3.1) 



as in the previous section. The correlation length i € of the boundary fluctuations is similar to tq, 

t £ =e- 1 i^ r0i (3.2) 

so that there is a non-trivial interaction between the boundary fluctuations and the wavefield. Here 
I is the scaled order-one correlation length defined below. 
2. The scale L £ of the propagation distance is much larger than Aq. More precisely, 



L e /A - 0{e~ 2 ). 



(3.3) 



Recall that the Rayleigh length for a beam with initial width r$ and central wavelength Ao is of the 
order of r^/Ao ~ £~ 2 Ao in absence of random fluctuations. The high-frequency scaling assumption 
(3.3) ensures that the propagation distance is similar to the Rayleigh length. 
3. The amplitude of the boundary fluctuations is small, of the order of £ 3 / 2 Ao. As we will show, this 
scaling is precisely the one that gives a cumulative scattering effect of order one after the propagation 
distance L £ . 

We use the hyperbolicity of the problem to truncate mathematically the boundary fluctuations to the 
range interval (Q,L/e 2 ). The bound L/e 2 is the maximum range of the fluctuations that can affect the 
waves up to the observation time T e of order e~ 2 . The lower bound in the range interval coincides with the 
location of the source. It is motivated by two facts: First, we observe the waves at positive ranges. Second, 
the backscattered field is negligible in the scaling regime defined above, as we show later in section |4~3} 

The boundary fluctuations are modeled with a random process ji 



T £ (x,z)=V l + e 3/2 fi(ex,ez) 



z e (0,L/s 2 ). 



(3.4) 
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The process /i is bounded, zero-mean, stationary and mixing, meaning in particular that its covariance is 
integrabl^] Because our method of solution flattens the boundary by changing coordinates, we require that 
/j, is twice differcntiable, with almost surely bounded derivatives. Its covariance function is given by 



and we denote by i? Q (£) its integral over £, 



Ro(0 



(3.5) 



(3.6) 



Our assumption on the differentiability of fi implies that R a is four times differentiable. Note that £ = is 
the maximum of the integrated covariance i? Q (£), so we have 



R' o (0) = 0. 



(3.7) 



We define the scaled square amplitude a 2 and correlation length I of the boundary fluctuations through the 
equations 



R o (0) = <J 2 



R>>(0) 



1 



(3.8) 



Ro(0) 

3.1. Change of coordinates. We introduce the change of coordinates from (x,y,z) to (x,rj,z), with 

yV 



n 



T £ (x,z)' 



(3.9) 



It straightens the boundary y — T £ (x,z) to r) = V, for any i£l and z £ (0,L/e 2 ). The pressure field in 
the new coordinates is denoted by 



P{u,x,T),z) = p\co,x, — ,z 

It satisfies the simple boundary conditions 

P{oj, x, T>, z) = 9^P(a;, x, 0, z) = 0, 

and the partial differential equation 



(3.10) 



(3.11) 



dl + dl 



V 2 

9 + V J 9 



d 2 - 2 V - 



2rj 



!VT 



e|2 



d v + k 2 



P 



f(oj,x, V )6(z), (3.12) 



derived from (1.1 1 and (3.10) using the chain rule. Here V and A are the gradient and Laplacian operators 
in (x, z) and f £ is the source of th e for m ( 2.15 ). 

When substituting the model (3.4) in (3.12), we obtain that P satisfies a randomly perturbed problem 

®x + 9* + (l - 2e 3/2 /x(ea;, ez)\ d 2 + k 2 + . . . P(w, x, rj, z) = f(u, ex, T))5(z). (3.13) 

The higher-order terms denoted by the dots are 

-2e 5/2 \l + 0(e 3/2 )] ??V/i • V0„P + 3e 3 fi 2 [l + 0(e 2 )] d 2 P - e 7/2 Afi \l + 0(e 3/2 )] rjd^P. 

They come from the expansions in e of the coefficients in (3.12), and are negligible in the limit e — > 
considered in section HI 



1 More precisely, /i is a ip-mixing process, with <p £ L 1//2 (R+), as stated in |13l 4.6.2]. 
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3.2. Wave decomposition. Equation (3.13) is not separable, but we can still write its solution in the 



L 2 (0,2?) basis of the eigenfunctions (2.7 1. The expansion is similar to (2.5) 

N 

P(w,x,T],z) = 'J2(t> j (r])u j (u,x,z) + 4>j(v)vj(u,x,z). 

3 = 1 3>N 

We define the forward and backward going wave mode amplitudes aj and bj by 

Oj(w,x,2) = (^-Uj(u,x,z) + 2i p ^ dzu 3 (u,x,zfj e~ l0 i ( " )z , 
bj(u),x,z) = (^Uj((j,x,z) - ^r^^ d z Uj(ui,x, zjje 1 ^^^ , 

so that the complex valued amplitudes of the propagating modes can be written as 

Uj(u, x, z) = aj {(J, x, z)e l ^ {u)z + bj(u, x, z)e~ ip ^ z . 



(3.14) 



(3.15) 



Definition (3.15) implies that 

d z a 3 (w, x, z)e if> ' {U1)Z + d s bj {u,x, z)e~ ip ' {ui)z =0, j = l,...,N. 



(3.16) 



This equation is needed to specify uniquely the propagating mode amplitudes, because they each satisfy a 
single boundary condition in the range (0,L/e 2 ) of the fluctuations. To derive these boundary conditions, 
let us observe that aj and bj must be constant in z G (—00, 0) and in z G {L/e 2 , 00), because the boundary 
is flat outside (0,L/e 2 ). Moreover, the radiation conditions 



lim dj(ui,x,z) = 0, lim bj(ui,x,z) = 0, 

z— > — 00 ' z— >oo 



imply that the mode amplitudes satisfy 



CLj{u), X, z = ) = 0, 

bj(uj, x, z = L/e 2 ) = 0. 



(3.17) 
(3.18) 



The last equation is the boundary condition for bj. The boundary value aj(ui,x, z = + ) follows from the 
jump conditions across the plane z = of the source in equation (3.13). We have 

[Uj]°- = 0, [d z Uj]°l = Fj{uj,ex), 



with Fj defined by (2.12). This gives 



[oj + bj]°- = 0, i/3j[aj - = F s {u, ex), 



and therefore 



aj(ui,x,0 + ) 



Fj{u>, ex). 



(3.19) 



Substituting (3.14) in (3.13), and using the orthonormality of the eigenfunctions <pj, we find that the 
wave mode amplitudes solve paraxial equations coupled by the random fluctuations in z G (0, L/e 2 ), 



(2ipjd z + d 2 x ) dj + t 2 ' ' '():>>, ra e 3/2 n{ex,ez)e- 1 ^ 



(-20jd z + flg) bj + e 2l ^ z d 2 x a 3 « e 3/2 ^{ex,ez)e l ^ z 



N 



qjivi 



U=l 

N 



1>N 



]T Iji (aie^' z + he~^ z ) + £ q,m 



2=1 



1>N 



, (3.20) 
(3.21) 



We dropped the higher-order terms that do not play a role in the limit e — > 0, and replaced the equality 
with the approximate sign. To simplify our notation, we omit henceforth all the arguments in the equations, 
except those of \l. The arguments will be spelled out only in definitions. 



The coupling matrix in (3.20 (-(3.21 ) is given by 



(3.22) 



It takes this simple diagonal form because we assumed a homogeneous background speed c . If we had a 
variable speed c(y), the matrix {qji} would not be diagonal, and the modes with j ^ I would be coupled. 
However, the results of the asymptotic analysis below would still hold, because the coupling would become 
negligible in the limit e — > considered in section [4j due to rapid phases arising in the right hand sides of 
d3720b, (13721 ). 



The equations for the evanescent components are obtained similarly, 

(d 2 + 8 2 - p]) Vj « e 3/2 fi(ex, ez) q^Vj, 



(3.23) 



and they are augmented with the decay conditions %(w, x, z) — > as s/x 2 + z 2 — > oo, for all j > N + 1, 

4. The limit process. We characterize next the wave field in the asymptotic limit e 0. We begin 
with the paraxial long range scaling that gives significant net scattering, and then take the limit. The scaling 
has already been described at the beginning of section [3] 

4.1. Asymptotic scaling. We obtain from ( 3.20[ )-( [3~22 ) that the propagating mode amplitudes satisfy 
the block diagonal system of partial differential equations 



20jd z 



e -2 lP] z d 2 

-2ip 3 d z + dl 



s 3 / 2 qjj[i(ex, ez) 



1 



e -2iP jZ 
1 



(4.1) 



for j — 1, . . . , N. Again, the appro xima te sign means equal to leading order. 

Because the right hand side in (4.1 ) is small, of order e 3 / 2 , and has zero statistical expectation, it follows 
from [5J Chapter 6] that there is no net scattering effect until we reach ranges of order e~ 2 . Thus, we let 



= Z/e 2 



(4.2) 



with scaled range Z independent of e. The source directivity in the range direction suggests observing the 
wavefield on a cross-range scale that is smaller than that in range. We choose it as 



X/e, 



(4.3) 



with scaled cross-range X independent of e, to balance the two terms in the paraxial operators in (4.1). 



regime (4.2)-(4.3). We denote them by 



Our goal is to characterize the e — > limit of the mode amplitudes in the paraxial long range scaling 
a E j(uj,X,Z) = aj (w,j, and b%uj,X,Z) = bj(u, ~ ~), (4.4) 



and obtain from (4.1 )-(|4.3[) that they satisfy the scaled system 



e -2iP 3 Z/s 2 Q2_ 

-2i^d z + d\ 



Z\ 




for j = 1, 



. , N, with initial conditions 
and end conditions 



2ipj(u) 



Fj(u,X), 



(4.6) 



b e Au),X,L) = 0. 



(4.7) 



4.2. The random propagator. Let us rewrite (4.5) in terms of the random propagator matrix 
P £ (w, X, X',Z)e C 2Nx2N , the solution of the initial value problem 



d z v e (u,x,x\z) = 



1 



m V [X,- H l u ,X, 



z 



G 



P e (w,X,X\Z), Z>0, 



P e (u,X,X',0) = 5(X-X')I. 



(4.8) 



Here I is the 2N x 27V identity matrix, S(X) is the Dirac delta distribution in X, and G and H are matrices 
with entries given by partial differential operators in X , with deterministic coefficients. We can define them 
from (4.5| once we note that the solution 

a\{u,X,Z) 



a e {uj,X, Z) = 



6 £ ( W ,X, Z) = 



a e N {uj,X, Z) 



t bf(u,X,Z) 
\ b%{uj,X,Z) 



(4.9) 



follows from 



a £ (w,X, Z) 
6 £ ( W ,X, Z) 



= j dX' P £ {lo,X,X' ,Z) 



a e (uj,X',0) 
b e {uj,X',Q) 



(4.10) 



Here b e (u>, X' ,0) is the vector of backward going amplitudes at the beginning of the randomly perturbed 
section of the waveguide, and it can be eliminated using the boundary identity 



a e (w,X,L) 




J dX' P £ (u},X,X',L) 



a s {u},X',0) 
b e (w,X',Q) 



(4.11) 



The initial conditions a e (w,X',Q) are given in (4.6). 



We obtain from (|4.5|) that H and G have the block form 

G = 



H 



H Q H b 

W< H 77 



G a Qb 

& G" 



where the bar denotes complex conjugation. The blocks are diagonal, with entries 

"V </..,- ,, , _ iSjiljj 



H 



2ft 



2ft 



and 



pa _ o2 pfc _ -2i0jZ/e 2 fP ■ j _ i 



Tjl 2ft Xl 11 2ft 



N. 



(4.12) 



(4.13) 



(4.14) 



The entries of the diagonal blocks depend only on the mode indices and the frequency, via ft(w). The entries 
of the off-diagonal blocks are rapidly oscillating, due to the large phases proportional to Z/e 2 . 

The symmetry relations satisfied by the blocks in H and G imply that the propagator has the form 



P £ (w,X,X',Z) 



T 6 (uj,X,X',Z) K^(uj,X,X',Z) 



R £ (w,X,X',Z) T e (w,X,X',Z) 
with N x N complex, diagonal blocks T £ and R £ . 



(4.15) 



4.3. The diffusion limit. The limit of P £ as e —> is a multi-dimensional Markov diffusion process, 
with entries satisfying a system of Ito-Schrodinger equations. This follows from the diffusion approximation 
theorem [14l[T5], see also [3 Chapter 6], applied to system (4. 



When computing the generator of the limit process, we obtain that due to the fast phases in the off- 
diagonal blocks of H and G, the forward and backward going amplitudes decouple as e — > 0. This implies 

9 



that there is no backscattered field in the limit, because the backward going amplitudes b £ are set to zero 
at Z = L. Equation (4.10) simplifies as 



a £ {uj,X, Z) = J dX'T £ (uj,X,X',Z)a £ (uj,X',0), 



where the initial conditions a £ (u), X' ,0) are given in (4.6). We call the complex diagonal matrix 
T £ (lo, X, X', Z) = diag (7?(w, X,X',Z),..., 7%(w, X, X', Z)) 



(4.16) 



(4.17) 



the transfer process, because it gives the amplitudes of the forward going modes at positive ranges Z ', in 
terms of the initial conditions at Z = 0. The limit transfer process is described in the next proposition. It 
follows straight from [T4l I15j . 

Proposition 4.1. As e — > 0, T £ (uj,X, X',Z) converges weakly and in distribution to the diffu- 
sion Markov process T(w, X, X', Z). This process is complex and diagonal matrix valued, with entries 
7j(oJ,X,X',Z) solving the Ito-Schrodinger equations 



dTj(uj,X,X\Z) = 



Q 2 _ 1 2 n R o{ Q ) 



2/?»^ 8(3] (u) 



Tj(u,X, X',Z)dZ 



■Tj{oj,X,X',Z)dB(X,Z), (4.18) 



for Z > 0, and initial conditions 



T j (u,X,X',0)=6(X-X'), j = l, 



Equations H-18) are uncoupled, but they are driven by the same Brownian field B(X, Z), satisfying 



E [B(X, Z)] =0, E [B(X, Z)B(X\ Z')\ = min{Z, Z'}R (X - X'), 



(4.19) 



(4.20) 



with R Q defined in (3.6). Thus, the transfer coefficients Tj are statistically correlated. 



The weak convergence in distribution means that we can calculate the limit e — > of statistical moments 
of T e , smoothed by integration over X against the initial conditions, using the Markov diffusion defined 



by (4.18 )-(4.19 ). In applications we have a fixed £ « 1, and we use Proposition 4.1 to approximate the 



statistical moments of the amplitudes of the forward going waveguide modes. 



When comparing the Ito-Schrodinger equations (4.18) to the deterministic Schrodinger equations (2.21) 



satisfied by the amplitudes in the ideal waveguides, we see that the random boundary scattering effect 
amounts to a net diffusion, as described by the last two terms in (4.18). We show next how this leads to 
loss of coherence of the waves, that is to exponential decay in range of the mean field. We also study the 
propagation of energy of the modes and quantify the decorrelation properties of the random fluctuations of 
their amplitudes. 



5. Statistics of the wave field. We begin in section [5TT] with the analysis of the coherent field. Ex- 
plicitly, we estimate the mean forward going mode amplitudes in the paraxial long range regime. Traditional 
imaging methods rely on these being large with respect to their random fluctuations. However, this is not the 
case, because E [of(o;, X, Zj\ decay exponentially with Z, at rates that increase monotonically with mode 
indices j. The second moments of the amplitudes do not decay, but there is decorrelation over the modes and 
the frequency and cross-range offsets, as shown in sections |5.3| and |5.4| Understanding these decorrelations 
is key to designing time reversal and imaging methods that are robust at low SNR. Robust means that wave 
focusing in time reversal or imaging is essentially independent of the realization of the random boundary 
fluctuations, it is statistically stable. Low SNR means that the coherent (mean) field, the "signal" , is faint 
with respect to its random fluctuations, the "noise" . 



5.1. The coherent field. The mean modal amplitudes arc 



e[o>,a:, z)} 



dX' E [Tj (w, X, X', Z)] a jM (w, X') , 
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(5.1) 



with mean transfer matrix satisfying the partial differential equation 



d z E[Tj(u,X,X',Z)] 



1 d 2 - 1 



with mode-dependent damping coefficients 



E[Tj{u,X,X',Z)], Z>0, 



8^H 



IV 2 \{N + a{u))-l/2) 2 -(J -f/2) 2 
(J-1/2) 4 



(5.2) 



(5.3) 



with units of length. Here we used definitions (2.13), (2.141 and (3.8), and obtained equation (5.3) by taking 

(5.4) 



expectations in (4.18). Its solution is given by 
E[Tj{u,X,X',Z)} 



2-rtiZ 



cxp 



Z 



2Z 



and the mean modal amplitudes are obtained from equations (5.1) and (4.6) 
E[a%u;,X,Z)] 



1 



2W 2n Pj(w)Z 
a,j t0 (lo, X, Z) exp 



dX' Fj(uj,X')exp 
Z 



Z 



t\2 



iP^jX- X') 
2Z 



(5.5) 



with a^ the solution of the paraxial wave equation ( 2.2l|2~22 ) in the ideal waveguide. 
The mean wave field follows from (3.14), after neglecting the evanescent part, 



E 



-/ X Z\~\ N 



Z -a t \ Z 



(5.6) 



It is different than the field in the ideal waveguides 



/ X Z\ N 
° ( w ' "7' r] ' ^2 J ra ^j(v) a ],o (w, X, Z) exp 



(5.7) 



because of the exponential decay of the mean mode amplitudes, on range scales Sj(ui). 

5.2. High-frequency and low-SNR regime. We call the length scales Sj(w) the mode- dependent 
scattering mean free paths, because they give the range over which the modes become essentially incoherent, 
with low SNR, 



SNRo 



\E[a%uj,X, Z)}\ 



E[\a%u,X,Z)\->] -\E[afa,X,Z) 



exp 



Z 



< 1, if Z^>Sj(u). 



(5.8) 



The second moments E |"|a f (a>, X, Z)| 2 ] are calculated in the next section, and they do not decay with range. 
This is why equation (5.8) holds. 

The scattering mean free paths decrease monotonically with mode indices j, as shown in (5.3). The first 
mode encounters less often the random boundary, and has the longest scattering mean free path 

SxM = A i^ p [{N + a(«) - 1/2) 2 - 1/4] » (5.9) 

The highest indexed mode scatters most frequently at the boundary, and its scattering mean free path 

2X> 2 a(uj) (2N + a(ui) - 1) a(u) S^uj) 



S n (uj) = 



a 2 ir 2 e (7V-1/2) 4 
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8 N 5 



(5.10) 



is much smaller than Si(oj), when N is large. To be complete, we also have 



Sj(uj) w5i(w)— j 



1-s 4 1 



iV 4 ' 



Hj=[aN\, ae(0,l), 



and 



(2J-1) 4 ' 



if j = o(iV). 



Our analysis of time reversal and imaging is carried in a high-frequency regime, with waveguide depth 
V much larger than the central wavelength A D or, equivalently, with N 1. We also assume a low-SNR 
regime, with scaled range Z exceeding the scattering mean free path of all the modes, so that none of the 
amplitudes <Zj are coherent. This is the most challenging case for sensor array imaging, because the wave 
field measured at the sensors is essentially just noise. We model the low-SNR regime using the dimensionless 
large parameter 



7 Si(wo) > ' 



(5.11) 



and observe from (5.3) that 



SjiuJo, 



>7>1, for all j = 1,...,N. 



(5.12) 



5.3. The second moments. The quantification of SNR and the analysis of time reversal and imaging 
involves the second moments of the mode amplitudes. Recall that 



a){u, X,Z)n J dX'Tf («, X, X\ Z) a jM (u, X') 



(5.13) 



with Tj the entries of the diagonal transfer matrix T £ . To calculate the second moments, we need to 

estimate E [TfTf] . The equations for Tf {ui x ,X X ,X[, Z)T t £ (uj 2 , X 2 , X^ Z) follow fr om the forward scattering 
approximation of (4.8 1, 



2&(wi) Al 2/3, (w 2 ) 



T 6 T £ 
'j 'i 



2eV2 



qjj fi(Xi,Z/e) qu fx(X 2 ,Z/e) 



A(W2) 



h 'i ' 



(5.14) 



for Z > 0, with initial condition 



rn^i,X 1 ,X' 1 ,0)Tl s {u} 2 ,X 2 ,X' 2 ,0) = - X[)5{X 2 - X' 2 ) 



(5.15) 



Their statistical distribution is characterized in the limit e — > by the diffusion approximation theorem 
OUS], see also [H Chapter 6]. It is the distribution of 7j(u>i, Xi, Z)Ti{w 2 , X 2 , Z), with Tj the limit transfer 
coefficients in Proposition |4.1| This gives the approximate relation 



E 



a e Au} 1 ,X 1 ,Z)af(ui 2 ,X 2 ,Z) 



xE 



dX[ J dX' 2 a j>ini (wi, X[) a tMi {u 2 ,X' 2 ) 
(ui , X x , X[ , Z)Ti(u 2 ,X 2 ,X 2 ,Z) 



(5.16) 



The calculation of E \j~jTi\ is given in appendix [a] We summarize the results in Propositions 5.l|5.3 
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5.3.1. The single mode and frequency moments. It is easier to calculate the diagonal moments, 
with j = I, and the same frequency ui\ = ui 2 = u). We have the following result proved in appendix [Aj 

Proposition 5.1. For all j = 1,...,N, and all the frequencies uj G [ujo — nB, uj + ttB], 



3 



T j (uj,X 1 ,X' 1 ,Z)T j (u;,X 2 ,X' 2 ,Z) 



2irZ 



exp 



t[3 J (u } )[(X 1 -X[f-(X 2 X^ 
2Z 



2Z 



■ f ds C [(Xi - X 2 )s + (X[ - X' 2 ){1 - s)] 
Jo 



with kernel C Q defined by 



C {X) = 1 - 



Ro(X) 
Ro(0) ' 



(5.17) 



(5.18) 



The general second moment formula does not have an explicit form in arbitrary regimes. But it can be 



approximated in the low-SNR regime (5.11). The expression (5.17) also simplifies in that regime, as stated 



in the following proposition, which we prove below. 



hand side in {5.17) is essentially zero, unless 

\X 1 -X 2 \ 



Proposition 5.2. In the low-SNR regime (5.11), and under the assumption X[ = X' 2 = X' , the right 



< 



7 Si (w) 



«1, 



and the moment formula simplifies to 



E 
with 



T j (u>,X 1 ,X',Z)T j (u>,X 2 ,X',Z) 



2itZ 



exp 



if3 j [(X 1 -X') 2 -(X 2 -X') 2 } {X x -X 2 ) 



X dJ (oj) = l\i 



2Z 



2XIM 



2Z ^y2 7 Si(w) 

If the initial points X[ and X 2 are different, but still close enough to satisfy 

\X[-X' 2 \ 



the moment formula becomes 



E 



T j {u,X 1 ,X[,Z)T 3 {u:,X 2 ,X 2 ,Z) 



2ttZ 



< 1, 



exp 



if3 j [{X 1 -X' 2 f-{X 2 -X' 2 f^ 



2Z 



x exp 



(X 1 - X 2 f + {X[ - X' 2 f + (X, - X 2 )(X[ - X' 2 ) 



Proof. We see from definitions (3.6 1 and (5.18) that C a [X) w 1 for |X| > t. Therefore 



/ ds C [(X x - X 2 )s] a 1 if \Xi - X 2 \ > £, 
Jo 



(5.19) 



(5.20) 



(5.21) 



(5.22) 



(5.23) 



and the right hand side in (5.17) becomes negligible, of order O (e 2Z / S i) <g; 1. In the case \X\ — X 2 \ ~ I 
we obtain similarly that the damping term is of order Z, and the right hand side in (5.17) is exponentially 
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small. It is only when \Xi — X 2 \ <C £ that the moment does not vanish. Then, we can approximate the 
kernel C in the integral with its first nonzero term in the Taylor expansion around zero, using the relations 



C o (0)=0, C£(0) = 0, andC;'(0) = - 



1 



i?o(0) 



(5.24) 



that follow from (|3.8[)-(|3.7|). We have 



n 



dsC [(X 1 -X 2 )s] 



|*i -X 2 I 2 



and the right hand side in 



5.17) is of the order exp 



U 2 
\x^x 2 \ 2 z 



if |X 1 -X 2 |»^, 



. This gives the condition (5.191, and the 



simpler moment formula (5.20) follows. 



Essentially the same proof applies in the case X[ ^ X 2 , because we can still expand the integrand in 



(5.17) by assumption (5.22). □ 

5.3.2. The two mode and frequency moments. The general second moment formula is derived in 
appendix [A] in the low-SNR regime (5.11 1. It has a complicated expression that we do not repeat here, but 
it simplifies for nearby frequencies, as stated below. 



Proposition 5.3. The modes decorrelate under the low-SNR assumption (5.11) 
e\t j (u 1 ,X 1i X[,Z)Ti(lo 2 ,X 2 ,X 2i Z)] «0 ifj^l, 



(5.25) 



for any two frequencies uj\ , w 2 and cross-ranges X\ , X 2 . The modes also decorrelate for frequency offsets that 
exceed 



S»/?f(u;)£ 2 _ p^u) ^-( w )/3»£ 2 



Z 2 \(3>{u)\ |/3»| 7 2S 2 H ' 
where f3j(uj) is the derivative of fij{oj) with respect to uj. For much smaller frequency offsets satisfying 

UJl + uj 2 



(5.26) 



\ui -u 2 \ < ^d,j (w), w 



(5.27) 



the moment formula is 



E 



Tj (oj 1 ,X 1 , X[, Z) Tj (wa, X 2 ,X' 2 , Z) 



/?» / i [p j {uj l ){X l - X[f - /3> 2 )(X 2 - X' 2 f 
exp < — 



2ttZ 



2Z 



{X 1 - X 2 f + (X[ - X' 2 f + (X 1 - X 2 ){X[ - X' 2 ) 
2X 2 H 



(5.28) 



5.4. Decorrelation properties. We already stated the decorrelation of the modes in Proposition [573] 
But even for a single mode, we have decorrelation over cross-range and frequency offsets. 



The decoherence length of mode j is denoted by Xd^(uj), and it is defined in (5.21 ). It is the length scale 



over which the second moment at frequency w decays with cross-range. It follows from (5.21) that Xjj is 



much smaller than the correlation length, for all the modes, and that it decreases monotonically with j. The 
first mode has the largest decorrelation length 



X d . 1 {uj)=t 



(5.29) 



because it scatters less often at the boundary. The decoherence length of the highest mode is much smaller 
in high-frequency regimes with TV 3> 1, 



Xd,N(u) — XdS 




N~ 5/2 . 



(5.30) 



The decorrelation frequency is derived in appendix 



A.2 



It is given by (5.26) or, more explicitly, by 



(N + a^-lf-ij-iy 



5/2 



64 7 2 2V" 9 (j - 1/2) 4 

it is much smaller than uj for all the modes, and it decreases monotonically with j, starting from 

wa 2 7r 3 0?/A) 3 



4 7 2jy4 



(5.31) 



(5.32) 



6. The forward model. Let us gather the results and summarize them in the following model of the 
pressure field 



/ X Z\ ^ 



dX'Tj(uj,X,X',Z) / drfhWffaX'rf), 



(6.1) 



where the symbol ~ stands for approximate, in distribution. That is to say, the statistical moments of 
the random pressure field P are approximately equal to those of the right handside. The first and second 
moments follow from Propositions |4.1j|5.3| In our analysis of time reversal and imaging we take small 
frequency offsets, satisfying |u5| <C £ld,j(w), so that we can use the simpler moment formula (5.28). 

The computation of the fourth moments of the transfer coefficients is quite involved. We estimate in 
appendix [B] some of them, for a particular combination of the mode indices and arguments. These moments 
are used in the next sections to show the statistical stability of the time reversal and coherent interferometric 
imaging functions. 

We analyze next time reversal and imaging in the low SNR regime, and assume for convenience that the 



source (2.15) has the separable form 



f(t,x, v ) 



Ox On 



X -X* 



1} — rj 



'v 



(6.2) 



meaning that the same pulse tp(t) is emitted from all the points in the support of the non- negative source 
density p. We scale this support with the dimensionless parameters Ox and V , and normalize the source by 



dX' 
~0x~ 



drj fX-X* 17 — 77* 



'A' 



= 1. 



The coefficients 



FAu,X) 



ixo v Jo 



V 



dr]4>j{ri)p 



X-X* 97 — 77* 



Ox 



(6.3) 



(6.4) 



are proportional to the Fourier coefficients (p(u)) of the pulse, and are thus supported in the frequency 
interval [w — B/2,u> + B/2]. The bandwidth B is small enough so that we can freeze the number of 



propagating modes to that at the central frequency, as explained in section 2.1 The width of the pulse ip(t) 
is inverse proportional to B, and we distinguish two regimes: The broadband regime with B 3> £ 2 w Q , and 
the narrowband regime with B < e 2 uj . The comparison with e 2 is because the source is at range Zj^/e 2 
from the array, and the modes arrive at time intervals of order 1/e 2 . Broadband pulses have smaller support 
than these travel times, meaning that we can observe the different arrivals of the modes, at least in the ideal 
waveguides. 

To analyze the resolution of time reversal and imaging, we study in detail the case of a source density 
localized around the point (X*, rf , 0). We say that we study the point spread time reversal and imaging 
functions, because the source has small support. Note however that it is not a point source. Its support is 
quantified by the positive parameters Ox and 0^ which are small, but independent of e. 
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7. Time reversal. Let us denote by D(t,X,T]) the pressure field measured in a time window ijj(t/T e ) 
at an array A, with aperture modeled by the indicator function 



l A (X,r 1 ) = l Ax (X)l Ari (r ) ), 



(7.1) 



at range z A = Z A /e 2 . Here X is the scaled cross-range in the array, related to the cross-range x by x — X/e, 
and Ax C K and A v C [0, T>] are intervals in X and 77. The window ip is a function of dimensionless 
arguments, of support of order one, and T £ denotes the length of time of the measurements. Because the 
waves travel distances of order e -2 , we scale T e as T £ = T/e 2 , with T of order one. 

In time reversal, the array takes the recorded field D(t, X, 77), time reverses it and emits D(T E — t, X, rf) 
back in the medium. We study in this section the resolution of the refocusing of the waves at the source, in 
the high-frequency and low-SNR regime described in section 5.2 Because we have a random waveguide, the 



resolution analysis includes that of statistical stability, given in section 7.3 



7.1. Mathematical model of time reversal. We have in our notation 



D(t,X,r)) = l A {X,r,)i>(~y(t, 



X 

e 



e 2 r 



(7.2) 



with mathematical model following from (6.1) 



D(t,X,rj)f*l A (X,ri)il> 
The time reversed field 

has Fourier transform 



3=1 



10 j (u 



2i(3j{uj) 



dX' F 3 {uj,X')T 2 {uj,X,X\Z a ). (7.3) 



D™(t,X,r J ) = D(T s -t,X, V ) 
D™(u>,X,r,) = e i " T 'D(u,X,r l ) 



(7.4) 
(7.5) 



with 



N 



D(cj,X,r,)^l A (X,r])J2<l>M 



i=i 



2n^ U ' 2iP 3 {uj - e 2 u/T) 



dX 



2 

£ U 



2 



,X')T J [u-— 1 X,X',Z A 



(7.6) 



The small frequency shifts e 2 u/T are due to the time scaling, and we can neglect them in the source terms 
Fj and in the amplitude factor 1//% . 

The model of the observed wave field at search locations (x s ,r/ s , z s ) — ( — , 77 s , is given by 



N 

o(t,x s , v s ,z s ) = Y / Mv s ) \ 

3=1 



dw 
27 



i/3j(ui) 



z A -z s 



iwt 



2ij3j(u)) 



x dX dr )( ( >j (ri)D' FR (uj,X,T,)T j (uj,X,X s ,Z A - Z s ), 



(7.7) 



using reciprocity. Note the similarity with equation (6.1), except that the source is now at the array, which 
we approximate in (7.7) as a continuum, instead of a discrete collection of sensors. This approximation 
is convenient for the analysis, because sums over the sensors are replaced by integrals over the X and r\ 
apertures, of lengths \Ax \ and \A n \. 
Using (7.5) in (7.7) and letting 



v 



dr)lAr,(v)<i>j{v)<l>l{'>l), 
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(7.; 



we obtain 

0(t,X', V ',Z') 



dX' f drf fX'-X* rf-rf 
p 



7X 



7 Y 



N 
3,1 = 1 



diu 
to* 



<p(u>) 



du ' 
2tt 



£(u) / dXl Ax (X)T J {uj,X,X s ,Z A -Z^Ti[ 



T ..V..V.Z.4 



x exp 



i/3j(w) 



e 2 u\Z A 



Tie 



(7.9) 



We define the time reversal function by 



j TE (r, jf) = o{t = t £ ,x s , v s 7 z s = o). 



(7.10) 



It models the wave field observed at the time instant t = T e , at the source range Z A . This is when and 
where the refocusing occurs. 

In the case of a source density that is tightly supported around (X*, rf), we may approximate J" rR by 



J TB (X s ,r, s ) 



duj 
2tt 



(7.11) 



with frequency-dependent kernel (point spread function) 



JY 



M™(u,X s ,r 1 s )^J2 T : 



3,1=1 



hWWiW) [du 



i/>(u) / dXl Ax (X) 



x exp 



Z A 



e 2 u\ Z A 



T 



(7.12) 



Here we used the source normalization (6.3). 

7.2. Resolution analysis. If the time reversal process is statistically stable, then we can estimate its 
refocusing resolution by studying the mean of (7.11). We refer to the next section for the analysis of the 
statistical stability of J"™. 

The mean time reversal function follows from (5.28) and ( 7.11 )-(7.12 ) 



E[J TR (X%ri s )} = J ^flu)E[M™(u,,X',if)], 



(7.13) 



with 



E[M TR (oj,X s ,r, s )} 



\Ax\ j^ Tjj ^ijfM^ ^(P'M Z A 



ldX 1 44^exp 



\A 



X - 



exp 
X s + X* 



{X s -X*) 2 
(X s -X*) 



(7.14) 



Moreover, letting 



A x = 



\Ax\ \A X \ 



2 
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(7.15) 



we obtain after integrating in X that 
E[M TR (oj,X s ,r, s )} 



8ttZ a 



- 33 



T 



exp 



*\2 



(X s - X*) 



X 



2Z A 



(X s - X*) 



exp 



2Z A 



[(* 



s\2 



*\21 



(X*) 



Note that 



(7.16) 



(7.17) 



are the scaled travel times of the modes, so only those modes that arrive within the support of the window 
ip contribute in (7.16|. 



7.2.1. Cross-range resolution. We observe in ( 7.16 I that modes contribute differently to the focusing 
in cross-range X, with resolution 



\X S -X*\< A XJ H := min X dd (w) 



2nZ A 



Recall from (5.21) and (5.29) that X^j decreases monotonically with j 

X dA (u3) [(N + a{w) 1/2) 2 - (j - 1/2) 2 ] 1/2 



whereas 



4(7 - 1/2)2 

2irZ A 
Pj{t»)\Ax\ 



N 



X dil (u)=l 



2Z A V 



[(N + a(uj)-l/2f-(j-l/2f]- 1/2 , 



(7.18) 



(7.19) 



(7.20) 



increases with j. Thus, in the high-frequency regime with N 3> 1, the cross-range resolution for the high- 
order modes is determined by the decorrelation length, even for large apertures. The cross-range resolution 
of the first modes may be determined by the aperture, but only if it is large enough, 



> 2Z a d J 2 1n. 



(7.21) 



It may appear at this point that the time reversal process can give good results even for small apertures 
\Ax\- However, we will see in section 7.3 that large apertures are needed for statistical stability. 

The modes with higher indices give the best cross-range resolution, but they travel at smaller speed. 
Thus, the focusing improves when we increase the recording time, because the array can capture the late 
arrivals of the high-order modes (see Figure 7.1). 

7.2.2. Depth resolution. To study the focusing in r), we evaluate the point spread function at cross- 
range X s = X* . We have 



3=1 



where Nt is the number of modes with arrival times in the recording window, 

Tj<T, for j = 1,2,..., N T < N. 

The coefficients Tjj are given by 



33 



v 



2( m \ _ \Ar)\ 



V 



m . 

— sine 
V 



2tt(7-1/2)| 
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m ■ 

— sine 
V 



27r(j-l/2)| 



>0, 



(7.22) 



(7.23) 



(7.24) 



Fig. 7.1. Depth profile (left) and cross range profile (right) of the mean point spread function for the time reversal 
functional. Here Z4 = 100, I = 1, a = 0.25, k = 60, T> = 1 (so that N = 19). The array diameter \Ax\ *s supposed to be 
smaller than the critical value \7.21]) which is about 220. Nt is the cut-off number (modes smaller than Nt are recorded and 
reemitted). Note that the high modes play an important role. The larger Nt is, the better the resolution. 



for an array in the set A v = [r]x, 772] C [0, T>]. They satisfy Tjj = 1 in the full aperture case A n — [0, T)\. 

The sum in (7.22 ) is maximum at rf = 77*, because all the terms are positive. The point spread function 
is smaller at other depths, because of cancellations in the sum of the oscillatory terms. We can make this 
more explicit in the high-frequency regime, with N ^> 1, if we write 
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and interpret (7.221 as a Riemann sum, which we then approximate with an integral. 
Consider for simplicity the full aperture case, where 
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The function A a becomes proportional to the Bessel function of first kind Jo as a 1, more explicitly, we 
have Ai(x) = (tt/2)Jq(x) so that 
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(7.27) 



We can then estimate the depth resolution as the distance between the peak of Jq, that occurs when r) s = if , 
and its first zero, that occurs when k\rj s —rj*\ w 2.4 (first zero of Jo). Therefore, the depth resolution of time 
reversal with full aperture is equal to the diffraction limit 



\r, s - V *\<A v (cj) 
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(7.28) 



if the array records the waves long enough to capture almost all the propagating modes. The resolution 
deteriorates if Nt is much smaller than N. Indeed for small a we have A Q (x) ~ (a/ 77 ) sinc(aa;) and therefore 
the depth resolution is A v (oj) « ttN/ (kNT) (see Figure 7.1). 
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7.3. Statistical stability. We now show that the time reversal function is statistically stable, meaning 
that the refocusing of the wave at the original source location does not depend on the the realization of 
the random medium but only on its statistical distribution, and the point spread function is approximately 
equal to its expectation. 

We restrict the analysis of statistical stability to the case of full aperture, where the calculations are 



simpler because the coupling matrix Tji becomes the identity. The point spread function follows from (7.12) 



and its variance at the source location is 

3? J — 1 

x{E[\T j (u,X,X*,ZX)\ 2 \Tj(u,Y,X*,Z A )\ 2 ] - E[\T 3 {u,X,X\Z a )\ 2 ]¥,[\Tj(lo,Y,X* ,Z A )\ 2 ]}. 

From appendix [B] (first case) we find that the variance is much smaller than the square expectation when 
\Ax\ £, and therefore the point spread function is equal to its mean approximately. The results contained 
in appendix [B] (second case) also show that if < ^ the variance of the point spread function is large, 
and therefore that time reversal refocusing may be unstable in this case. 

There is however another mechanism that can ensure statistical stability of the focal spot if the array is 
small. Indeed, if the bandwidth of ip is larger than the decorrelation frequency, then the variance 

V^[J^(X*,r,*)] = J^-J ^W^^Co^M^X^^M^'^X^rf)] 

is small because the covariance of the point spread function at two frequencies becomes approximately zero 
if the frequency gap is large enough. Therefore, if the pulse has large bandwidth, then the time-reversal 
focal spot is statistically stable even for small arrays. 

8. Imaging. The sharp and stable focusing of the time reversal process in the random waveguide is 
due the backpropagation of the time reversed field D TR in exactly the same waveguide. Time reversal is 
a physical experiment, where the waves can be observed in the vicinity of the source, as they refocus. In 
imaging we only have access to the data measured at the array, and the backpropagation to the search 
points is synthetic. Because we cannot know the fluctuations of the boundary, we simply ignore them in the 
synthetic backpropagation and obtain the so-called reversed time migration imaging function. We analyze 
it in section [O] and show that it does not give useful results in the low-SNR regime. In particular, we show 
that the images are not statistically stable with respect to realizations of the fluctuations. Stability can be 
achieved by imaging with local cross-correlations of the array measurements. Local means that we recall 
the decorrelation properties of the random mode amplitudes described in section [5^4} and cross-correlate the 
measurements over receivers located at nearby cross-ranges X, and projected on the same eigenfunctions. 
The resulting coherent inter jerometric imaging method is analyzed in section |8.2| 



8.1. Reverse time migration. The reverse time migration function is given by the time reversed 
data D TR propagated (migrated) in the ideal waveguide to the search points (x s , ij s , z s ) — (^-,7? s , fr) • Its 



mathematical expression follows from (2.20), with amplitudes (2.18) replaced by 



a jt0 (u,X,Z)~* I dX% (u,X,X',Z)--? T ^ [ dvhWD™^^',?]). (8.1) 



The ideal transfer coefficients 
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are defined by the Green's functions of the paraxial operator in (2.21 ). We obtain 
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with the right hand side evaluated at the same time t — T e as in time reversal. 

We as sume again a tightly supported source density normalized by (6.3) and substitute the model (7.5) 
of D TR in |8!3| to obtain 
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with frequency-dependent kernel (point spread function) 
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(8.5) 
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8.1.1. The mean imaging function. Let us take for simplicity the case of full aperture in depth, 
where the coupling matrix Tji given by (7.8) becomes the identity. We obtain from (8.5) and the moment 
formula (5.4) that 
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(8.6) 



Moreover, assuming the aperture Ax defined in (7.15), and integrating in X, we get 
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(8.7) 



This result is almost the same as in the ideal waveguide, except for the damping coefficients exp [— Z A /Sj] . 

The sine kernel in the mean point spread function gives the focusing in cross-range, with mode-dependent 
resolution 
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The best resolution is for the first mode, that has the largest wavenumber Px{uS) ps ttN/T> sa k, and gives 
the Raylcigh cross-range resolution 
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The focusing of the point spread function M M in range can only be due to the summation of the rapidly 
oscillating terms exp [— i(3jZ s /e 2 ] . But these terms are weighted by exp[—Z A /Sj], which decay fast in j. 
The first term dominates in 
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(8.10) 



so the mode diversity does not lead to focusing in range, as is the case in ideal waveguides. Nevertheless, 
the mean reverse time migration function peaks at Z s = because of the integral over the bandwidth in 



E[J M (X\ V \Z S )} = J ^(co)E[M M (u,X\ V \Z s )} 



(8.11) 



and the range resolution is of the order e 2 /[f3[(ui )B]. 

When we evaluate the point spread function at Z s = and X s = X* , we obtain 



E [AT>, X* = X\ if, Z-0)] = Mf tMMfl ^ ( e -*x/«,(«>. (8 . 12) 



This is a sum of the oscillatory functions 
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multiplied by positive weights, which are small and decay fast in j. The first term dominates in (8.12) and 
there is no depth resolution at all. We show next that these small weights also indicate the lack of statistical 
stability of the reverse time migration function. 

8.1.2. Stability analysis. To assess the stability of the reverse time migration, we calculate its variance 
at the source location 
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We have from the results above that 
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The second moment of J m is 
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and we recall from Proposition |5.3| that only the diagonal terms j = I contribute to the expectation. We also 
assume a small bandwidth B <C Qdj, for all the modes j, so that we can use the simpler moment formula 
|K28]. We obtain 
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after approximating the modal wavenumbers by their value at the central frequency. This expression can be 
approximated further, after integrating in X\ and X 2 , and supposing that the decoherence lengths Xdj are 
much smaller than the array aperture, 
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(8.16) 



The second moment (8.16) is clearly much larger than the square of the mean (8.13), which is expo- 
nentially small in range. Although the mean of the imaging function is focused at the source, it cannot be 
observed because it is dominated by its random fluctuations. The reverse time migration lacks statistical 
stability with respect to the realizations of the random fluctuations of the boundary of the waveguide. 

The calculations above are for a small bandwidth, satisfying B <C Cldj for all the modes captured in 
the recording window. The calculations are more complicated for a larger bandwidth, but the conclusion 
remains that reverse time migration is not stable with respect to different realization of the random boundary 
fluctuations. 

8.2. Coherent interferometric imaging. The main idea of the coherent interferometric (CINT) 
imaging approach is to backpropagate synthetically to the imaging points the local cross-correlations of 
the array measurements, instead of the measurements themselves. By local we mean that because of the 
statistical decorrelation properties of the random mode amplitudes described in section [5T4} we cross-correlate 
the data D(uj,X, 77) at nearby frequencies and cross-ranges X, after projecting it on the subspace of one 
cigcnfunction <pj at a time. The projection gives the coefficients 
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which are directly proportional to the coefficients Fj of the source only in the case of an array spanning the 
entire depth of the waveguide. We assume this case here, because it simplifies the analysis of the focusing 
and stability of the CINT function. We also take a small source, meaning that we essentially compute the 
CINT point spread function. 

The model of the coefficients ( |8.17 1 is 
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and we cross-correlate them at cross-ranges satisfying \Xi — X 2 \ < Xd.j(Lo), and at frequency offsets 
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We take such small Cl to simplify the second moment formulas. 

The CINT image is formed by backpropagating the cross-correlations to the imaging point, using the 
Green's function in the ideal waveguide. We first define the CINT image in the (X, Z)-domain: 
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(8.21) 



where lx d j are indicator functions of the cross-range interval [— Xdj(co), Xdj(uj)] calculated at the central 
frequency u> = {lo\ + ui 2 )/2. Similarly, 1q is the indicator function of the frequency interval [— CI, CI]. 
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8.2.1. The mean CINT function. To study the focusing of CINT, we consider its expectation 
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This expression follows from (8.21), the second moment formula (5.28), and definition (8.2) of the ideal 
transfer coefficients 7},o- 



8.2.2. Cross-range focusing. Let us consider in (8.23) a search point at the range of the source 
Z s = 0, 
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This formula simplifies after integrating over the array aperture and assuming as before that X^j <C \Ax 
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Each term in the sum focuses at the source, with resolution 
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defined as twice the standard deviation of the Gaussian in (8.25). The number of modes participating in the 
sum is determined by the length of the recording time window, as before, but each mode is weighted by the 
correlation length X&j, which decreases monotonically with j. The first mode has the largest contribution in 



5.25), and gives the best cross-range resolution. Since its wavenumber is approximately fa{oS) ~ irN/T> w k, 
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is comparable to the classic Rayleigh resolution for an array of aperture equal to ivXd,i(uj) (see Figure 8.1 ). 

The cross-range resolution (8.27) is worse than that of time reversal. Scattering at the random boundary 
is beneficial to the time reversal process, and the more modes are recorded, the better the result. However, 
scattering impedes imaging, and the best cross-range resolution is achieved with the first mode. Even with 
this mode, the resolution is worse than that in ideal waveguides 2irZ A / (k\Ax\), because X^ t \ <C |-4-x|- 
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8.2.3. Range focusing. When we evaluate the mean CINT point spread function (8.23) at the cross- 
range X s — X*, we obtain 
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Because we integrate over lj the rapidly oscillating integrand, at scale e 2 , we have from the method of 
stationary phase that (8.28) is large for 



Z s = s 2 C 



with ( s independent of e. Recall the assumption (8.19) of the frequency offsets. 
The mean point spread function becomes 



and we define the mode-dependent scaled range resolution by 
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Again, the resolution is best for the first mode, which has the largest weight X^i{u)) in (8.29). See Figure 
18.11 for an illustration. 



8.2.4. Depth estimation. One natural way to estimate the depth rf would be to consider the full 
CINT imaging functional 
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This is a sum of positive terms and it does not have a peak at the depth of the source (see Figure 8.2 ). 

Because of scattering at the random boundary the modes are decoupled, and we cannot speak of coherent 
imaging in depth. We work instead with the squares of the mode amplitudes, i.e. intensities. Incoherent 
imaging means estimating the depth of the source based on the mathematical model (8.32). More explicitly, 
we can estimate rj* by solving the least squares minimization problem 
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where the estimators X* and Z* of the cross-range X* and range offset Z* = of the source have been 
determined as the location of the maximum of (8.20) (see Figure 8.2). 

25 



Fig. 8.1. Range profile and cross range profile of the mean point spread function for the CINT functional. Here Zj± = 100, 
I = 1, a = 0.25, k = 60, T> = 1 (so that N = 19,), and the cut-off freguency is f!/c = 1. Nt is the cut-off number (modes 
smaller than Nx are recorded and reemitted). Note that the high modes do not play any role. 




Fig. 8.2. Depth profile with \8.31\ (left) and with ^8. 33\ l (right) for the CINT functional. In the right picture we plot the 
reciprocal of the square root of the function in \8.33\ . Here = 100, £ = 1, a = 0.25, k = 60, T> = 1 (so that N = 19). Nj< 
is the cut-off number (modes smaller than Nt are recorded and reemitted). Note that the high modes do not play any role. 



8.2.5. Statistical stability. The analysis of statistical stability of the CINT function is basically the 
same as that of time reversal. The function is stable when evaluated in the vicinity of the source location 



if the array has large aperture \Ax\ ^> We have seen in section 8.2.2 that a large aperture does not 
improve the focusing of E [J' ™ 17 ] . The cross-range resolution is limited by the decoherence length. But a 
large aperture is needed for the CINT function to be statistically stable. 

Another way of achieving statistical stability of CINT is to have a pulse with large bandwidth. This was 



already noted in the discussion of statistical stability of time reversal in section 7.3 



Note that the statistical stability of CINT relies on computing correctly the local cross-correlations of the 
measurements at the array. By this we mean that the cross-range and frequency offsets in the correlations 
should not exceed the decoherence length and frequency. Moreover, the cross-correlations should be with 
one mode at a time. This can be done with arrays that span the whole depth of the waveguide, because the 
coupling matrix B^ becomes the identity when \ A^\ — T>. If the aperture \ A„\ is small, there arc large mode 
index offsets \j — l\ for which Tji ^ 0. Consequently, there are many terms of the form 7} 71, with j =/= I, 
that participate in the expression of the imaging function. Since only the diagonal terms are correlated, we 
obtain that J CINT has large variance when \A V \ <§; V. 

In practice, the decoherence scales Xdj and fldj are likely not known explicitly. The formulas derived 
above are specific to our mathematical model. However, the decoherence scales can be estimated as we form 
the image, using an adaptive procedure similar to that introduced in 0. 
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9. Summary. In this paper we analyze propagation of acoustic waves in three-dimensional random 
waveguides. The waves are trapped by top and bottom boundaries, but the medium is unbounded in the 
remaining two directions. The top boundary has small, random fluctuations. We consider a source that 
emits a beam, and study the resulting random wave field in the waveguide. 

The analysis is in a long range, paraxial scaling regime modeled with a small parameter e. It is defined 
as the ratio of the central wavelength A D of the pulse emitted from the source and the emitted beam width 
r . The range of propagation is of the order of the Rayleigh length Tq/X = e~ 2 \ . The fluctuations of the 
boundary are on a length scale that is similar to the beam width, and their small amplitude is scaled so 
that they cause significant cumulative scattering effects when the waves travel at ranges of the order of the 
Rayleigh length. 

The wave field is given by a superposition of waveguide modes with random amplitudes. The modes are 
solutions of the wave equation in the ideal waveguide, with flat boundary. The scattering effects are captured 
by their random amplitudes. We show that in our scaling regime the amplitudes satisfy a system of paraxial 
equations driven by the same Brownian motion field. We use the system to calculate three important mode- 
dependent scales that quantify the net scattering effects in the waveguide, and play a key role in applications 
such as imaging and time reversal. The first mode-dependent scale is the scattering mean free path. It gives 
the range over which the mode loses its coherence, meaning that the expectation of its random amplitude is 
smaller than its fluctuations. The other mode-dependent scales are the decoherence length and frequency. 
They give the cross-range scale and frequency offsets over which the mode amplitudes become statistically 
uncorrelated. 

We use the results of the analysis to study time reversal and imaging of the source with a remote array of 
sensors, in a low SNR regime. Low SNR means that the waves travel over distances that exceed the scattering 
mean free paths of all the modes, so that the random wave field measured at the array is dominated by its 
fluctuations. 

In time reversal, the waves received at the array are time reversed and then re-emitted in the medium. 
They travel back to the source and refocus. The refocusing is expected by the time reversibility of the wave 
equation, but the resolution is limited in ideal waveguides by the aperture of the array. We analyze the time 
reversal process in the random waveguide and show that super-resolution occurs, meaning that scattering at 
the random boundary improves the refocusing resolution. An essential part of the resolution analysis is the 
assessment of statistical stability with respect to different realizations of the random boundary fluctuations. 
We show that statistical stability holds if the array has large aperture and/or the emitted pulse from the 
source has a large bandwidth. 

Time reversal is very different from imaging. In time reversal the array measurements are backpropagatcd 
physically, in the real waveguide. In imaging we can only backpropagate the time reversed data in software, 
in a surrogate waveguide. Because we cannot know the boundary fluctuations, we neglect them altogether, 
and the surrogate is the ideal waveguide. The resulting imaging function is called reverse time migration 
and it does not work in low SNR regimes. It lacks statistical stability, i.e., the images change unpredictably 
from one realization of the fluctuations to another. 

We show that robust imaging can be carried out in low SNR regimes if we backpropagate local cross- 
correlations of the array measurements, instead of the measurements themselves. Here local means that 
we cross-correlate the data projected on one mode at a time, and for nearby cross-ranges and frequencies. 
The method is called coherent interferometric (CINT), because it is an extension of the CINT approach 
introduced and analyzed in [5J 03 21 [2] for imaging in open, random environments. We show that CINT 
images are statistically stable under two conditions: The first condition is the same as in time reversal and 
it says that the array should have a large aperture and/or the pulse bandwidth should be large. The second 
condition is that the cross-range and frequency offsets used in the calculation of the local cross-correlations 
do not exceed the mode-dependent decoherence length and frequency, respectively. We derived mathematical 
expressions of these scales, for our model. In practice, they can be estimated adaptively, using the image 
formation, with an approach similar to that in [3]. The estimation is possible because there is a trade-off 
between stability and resolution that is quantified by the decoherence scales. If we over-estimate them we 
lose statistical stability. If we undcr-estimatc them, we lose resolution. 

While cumulative scattering aids in time reversal, it impedes imaging. We quantify this explicitly in the 
resolution analysis of CINT. In time reversal the resolution improves when we record the wave field over a 
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long time, so that wc include the high-order modes that travel at slower speed. In CINT, the best cross-range 
and range resolution is given by the first mode, which encounters the random boundary less often, and is 
thus less affected by the fluctuations. The cross-range resolution is similar to the classic Rayleigh one of 
range times wavelength divided by the aperture, but instead of the real aperture we have the decoherence 
length of the mode. This length decreases monotonically with range, because longer distances of propagation 
in the random waveguide mean stronger scattering effects. Similarly, the range resolution is similar to the 
classic one, of speed divided by the bandwidth, but the bandwidth is replaced by the decoherence frequency 
which decreases monotonically with range. 

The estimation of the depth of the source is different than that of range and cross-range. Because the 
modes decorrelate in the low SNR regime, we cross-correlate the data projected on one mode at a time, so 
essentially, we work with intensities. The estimation of the depth of the source from the intensities can be 
done by minimizing the misfit between the processed measurements and the mathematical model. While 
the cross-range and range estimation with CINT is done best with the first waveguide mode, the depth 
estimation requires many modes. Thus, we still need a long recording time at the array to capture the later 
arrival of the high-order modes. 
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Appendix A. Second moment calculation. 

The equation for E [Tj(ui ) Xi,X{,Z)Ti{u)2,X 2 ,X! 2) Z)\ follows from (4.18), using Ito calculus, 



d z E [73-77] 



-d 2 



X, 



2/3; M 



x. 



_ 2C (X 1 - X 2 ) 



E [757T] .(A.l) 



Its solution can be written as 



E 



T j (u)i,X 1 ,X' 1 ,Z)Ti(co 2 ,X 2 ,X 2 ,Z) =M j i(u 1 ,u 2 ,X 1 ,X 2 ,Z;X' 1 ,X! l )e 



with Mji solving 



d z Mji = 



i o2 i pp. 2C (A! - X 2 ) 



for Z > 0, and the initial condition 

M^uuut, X u A 2 ,0; X[, X' 2 ) = 5{X l - X[)S(X 2 - X' 2 ) 



Mi, 



(A.2) 



(A.3) 



(A.l) 



A.l. Single frequency. Let us begin with the single frequency case, Wi = uj 2 — w, and introduce the 
center and difference coordinates £ and £ so that 



In this coordinate system we have that 

U jl (u } ,t,lZ;e,e) = M jl (< 



U),U), 



.Z 



(A.5) 



(A.6) 



satisfies the initial value problem 



dzUji = nki),J 



--C 



1 1 
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1 1 A £ 



u jh Z>0, 



(A.7) 



Its Fourier transform in £ is the Wigner distribution 

the solution of the transport equation 

[d z +W t ]W jl {^,K,Z;£, 1 ,1') = 



2tt 



4VA 



x exp 



dqC U 



ii 



for Z > 0, with initial condition 



and kernel 



Here 



C (k) = (S(k) 



A. 1.1. Single mode moments. The transport equation | |A.9 | simplifies in the case j = I, 



[d z + W 6 }W jj (u J ,Z,K ) Z;e ) e) = - j- / dqC {q)W jl (oJ > ^K 



^■(w,«,e,z ; e',e') = 

satisfies the initial value problem 



W jj (cj,t,K,0;e,e) = ^e- i *5(Z-a 

in k and f. 



and can be integrated easily after Fourier transforming in k and £. Explicitly, 
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dz + KC?; 



^i(w,K,€,^;^,C ? ) = -| : C (-^=)K iJ (w,K,f,Z;C,C ? ), z - 



2ir 

which can be solved with the method of characteristics. 
We obtain that 



2tt 



kZ ) exp 



Sj Jo 



and tracing back out transformations (A.2), (A.6), (A.8) and (A.13), we get 

~ X' x ) 2 — (X 2 — X' 2 ) 2 ] 
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T j (uj,X 1 ,X' 1 ,Z)T j (u,X 2 ,X 2 ,Z) 



2nZ 



exp 
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This is the result stated in Proposition |5.1 



ds C a 
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A. 1.2. Two mode moments. It is not possible to obtain a closed form solution of (A. 9 1, unless we 
make further assumptions. We consider the low-SNR regime described in section 5.2 and suppose that 

\X!-X 3 \ <X dJ (u;)«£ (A.17) 

This is the condition under which the diagonal moments E [7} 77] are not exponentially small, by Proposition 



5.2 The two mode moments cannot be larger than the diagonal ones, so they are essentially zero when ( A.17 ) 



does not hold. 



Note that in (A. 9 1 n is the dual variable to £ ~ W~f3j(Xi — X 2 ), and that q is in the support of C Q , so 



\q\ < l/£. Therefore, 



%\x 1 -x 2 JWji VP 



and we can expand the Wigner transform in (A.9) around n. The exponential can also be expanded when 

(A.18) 



meaning that j and I are close. We return at the end of this section to this point 



Assumptions ( A.17 1-( A.18) justify the approximation of the right hand side in (A.9) by the second-order 



expansion in q of the product of the exponential and the Wigner transform. We obtain that 
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for Z > 0, with initial condition (A. 10). This equation is solved in [TJ. The result follows from the inverse 
Fourier transform in k of the solution, and from (A. 2), (A. 6), 
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Formula (A.20) is complicated, but it can be simplified under the assumption that 



Z 2 \P 3 -Pi\ 



H j M 2 ^/s~Si 



< 1. 



(A.20) 



(A.21) 



Then, we can expand the sine, cot and sin functions in (|A.20) and obtain the simpler formula 

~ i{Pj{Xi — X^) 2 — (3i(X 2 — X' 2 ) 2 )~ 
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It remains to justify assumptions (A. 18) and (A. 21 1. Because of the exponential decay in Z, we note 



that the moments are essentially zero unless 



Z< 1. 



But in our low-SNR regime this translates to 

('- 



< « 1, 
7<Si 



by definition (5.11), and it is satisfied only when j = I. This justifies the assumptions, and it means that 
the modes are essentially decorrelated. 

A. 2. Two frequency moments. The calculation of the two frequency moments is exactly as in the 
previous section, with j3j replaced by (3j(uii) and $ replaced by f3j(w 2 ). We only consider the case j = Z, 
because the modes decorelate as explained above. The moment formula follows from (A. 20), with j3j replaced 
by /3j(wi) and /3; replaced by (3j(uj 2 ), and similar for Sj and Si. We can simplify it under the assumption 
that \wi — ui 2 \ is sufficiently small to make first-order expansions in u)\ — u) 2 . Let uj and uj be the center and 
difference frequencies 



uj = u)\ - uj 2 , 



We have from \k.2Q\ that 
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(A.23) 



A. 3. Frequency decorrelation. To study the decorrelation over frequency offsets, let Xj — X 2 and 



X[ = X' 2 in (A.23) 
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We have two factors that decay exponentially in uj. The first is the sine, decaying at the rate 



and the second is the Gaussian with standard deviation 



\/3'An)\ 7 2 Si 2 H 



>2Z 
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\p'Au>)\ V ZyS^u) 



(A.24) 
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(A.26) 



Note that 
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(A.27) 



and using equations (2.13), (2.14), (5.3), and the high-frequency assumption N ^> 1, we have 
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and £/A = O(l), because the scaled correlation length is similar to the wavelength A. Moreover, the rate 
( A.26| is given by 



(iV + aM-I) 2 -(i-|)' 



3/2 
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and it is larger than Cl^ : j(u>). 

Thus, we call Q,dj(ui) the mode-dependent decoherence frequency, the frequency scale over which the 
second moments decay. Note that when the frequency offsets satisfy \uj\ -C &d,j the moment formula (A. 24) 
simplifies to expression (5.28) in Proposition 5.3 because 



exp 



2ft 2 (w) 
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when 



\u>\ < Q, dtj < sij(ui). 

Appendix B. The fourth moments. We denote the moments by 
M. jLJ l ■= E 
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and obtain from (|4.18|) that they satisfy the partial differential equation 
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for Z > 0, with the initial condition 

Mjiji, = 5{Xx - X[)S(X 2 - X' 2 )5{Y X - Y{)S{Y 2 - Y 2 '), at Z = 0. 
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Let us consider the case j = I, J = L, w\ = uj 2 = u>, and uj^ — oj^ = a/. These moments Mjjjj are 
needed in section [7] to show the statistical stability of the time reversal function in the case of an array that 
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spans the entire depth of the waveguide. We look for the fourth-order moment for X\ = X 2 and Y\ = Y 2 in 
the support of the array. So we parameterize 
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Equation (B.2) becomes (remember C"(0) = l/i 2 ) 
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with the initial condition 
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We address two cases: 

Case 1: The array diameter \Ax\ is much larger than t. This allows us to simplify Equation (B.8 1 as 
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which has a separable form in and (w,C)j ancl we § e t (following the same method as in the case of 

second-order moments): 
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Equivalently, in terms of the original variables, 
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which is equal to E [7^ 77] E [TjTj] ■ 

Case 2: The array diameter \ Ax\ is smaller than I, Then Equation (B.2) becomes 
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with the initial condition (B.9). This equation can be solved explicitly after Fourier transforming in £ and 
C If we let 
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The solution is given by the method of characteristics 
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and the moment estimate follows from the inverse Fourier transform, 
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Equivalently, in terms of the original variables, 
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If Xi = X 2 and Yi=Y 2 , then 
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Here we can see that the fourth-order moment is not equal to the product of the second-order moments. 
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